#pragma once
#include <vector>
using namespace std;
double divided_diff(vector<double> vx, double (*f)(double)) // compute the f[x0, x1, ..., xk]
{
	int k = vx.size() - 1; // this "k" is denoted the function will return a k-th divided difference 
	
	if (k == 0)
	{
		return f(vx[0]);
	}
	else
	{
		vector<double> vxtemp1 = vx, vxtemp2 = vx; //this design is a little bit stupid...
		//it1 --- f[x1, ..., xk]; it2 --- f[x0, ..., xk-1]
		vector<double>::iterator it1 = vxtemp1.begin(), it2 = vxtemp2.end() - 1;
		double x0 = *it1, xk = *it2;
		vxtemp1.erase(it1);
		vxtemp2.erase(it2);
		return (divided_diff(vxtemp1, f) - divided_diff(vxtemp2, f)) / (xk - x0);
	}
}
double PI_n(vector<double> vx, double x) // pai_n(x)
{
	int n = vx.size();
	if (n == 0)
	{
		return 1.0;
	}
	else
	{
		double result = 1.0;
		for (int i = 0; i < n; i++)
		{
			result *= x - vx[i];
		}
		return result;
	}
}
double newtonIP(vector<double> vx, double (*f)(double), double x)
{
	int n = vx.size() - 1;
	vector<double> a(n+1);
	for (int k = 0; k <= n; k++) 
	{
		vector<double> vxtemp1 = vx;
		vxtemp1.erase(vxtemp1.begin() + k + 1, vxtemp1.end()); // reserve from x0 to xk 
		a[k] = divided_diff(vxtemp1, f);
	}
	double result = 0.0;
	for (int k = 0; k <= n; k++)
	{
		vector<double> vxtemp2 = vx;
		vxtemp2.erase(vxtemp2.begin() + k, vxtemp2.end()); // reserve from x0 to x(k-1) 
		result += a[k] * PI_n(vxtemp2, x);
	}


	return result;
}
